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s Abstract 

9 Systems where resource availability approaches a critical threshold are com- 

10 mon to many engineering and scientific applications and often necessitate 

11 the estimation of first passage time statistics of a Brownian motion (Bm) 

12 driven by time-dependent drift and diffusion coefficients. Modeling such sys- 

13 terns requires solving the associated Fokker-Planck equation subject to an 

14 absorbing barrier. Transitional probabilities are derived via the method of 

15 images, whose applicability to time dependent problems is shown to be lim- 

16 ited to state-independent drift and diffusion coefficients that only depend on 

17 time and are proportional to each other. First passage time statistics, such 
is as the survival probabilities and first passage time densities are obtained 

19 analytically. The analysis includes the study of different functional forms of 

20 the time dependent drift and diffusion, including power-law time dependence 

21 and different periodic drivers. As a case study of these theoretical results, 

22 a stochastic model of water resources availability in snowmelt dominated re- 

23 gions is presented, where both temperature effects and snow-precipitation 

24 input are incorporated. 
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27 1. Introduction 



A wide range of geophysical and environmental processes occur under 
the influence of an external time-dependent and random forcing. Climate- 



30 driven phenomena, such as plant productivity (Ehleringer et al. , 1997), steno- 



thermal populations dynamics (McClanahan & Maina, 2003), crop produc 



32 tion (Rosenzweig & Parry, 1994), the alternation between snow-storage and 



33 melting in mountain regions (Hamlet & Lettenmaier, 1999 Marks et al. 



,1998 ), the life cycle of tidal communities (Barranguet et al. , 1998 Bertness 



& Leonard, 1997 Charles & Dukes, 2009), and water-borne diseases out 



breaks (Pascual et al. , 2002; Patz et al. 2005) offer a few such examples. In 



particular, several environmental systems can be described by state variables 
representing the availability of a resource whose dynamics is forced by diverse 
environmental factors and climatic oscillations. Elevated regions water avail- 
ability - mainly originating from the melting of snow masses accumulated 
during the winter period (under the forcing of increasing temperatures), and 
precipitation (moving from the solid precipitation to the rainfall regime) - 
offers a relevant case study (presented in Section |ij). All of these processes 
are now receiving increased attention in several branches of ecology, climate 
sciences and hydrology, due to their inherent sensitivity to climatic variabil- 
ity. 

Analogous dynamical patterns can be found in slowly-driven, non-equili- 
brium systems with self organized criticality (SOC), where the density of 



2 



potentially relaxable sites in the system can be described via a random 
walk with time-dependent drift and diffusion terms (Adami, 1995 Bak & 



Paczuski , 


1995 


Jensen 


1998) 



diffusion term derives from a gradual decrease of susceptible sites, so that 
sites availability acts on the directionality and pathways (drift term) of the 
"avalanches" till diffusion "kills" all the activity in the system ( Redner| |2001 



pp. 120-131). Similar dynamics occur in systems displaying stochastic reso- 
nance, where noise becomes modulated by an external periodic forcing (see 
Bulsara et aT| [19941 [T9951 |Gammaitoni et al.[ pTOj |McDonnell et al.[ |2008| 



and references therein). 

In many instances, the above mentioned processes are restricted to the 
positive semi-plane or to the time at which a certain critical threshold is 
reached, and are represented by a Fokker-Planck (FP) equation with an ab- 
sorbing barrier. The main focus here is on the first passage time statistics 
of the process, such as the survival probabilities and the first passage time 
densities. In the following, a brief review of the general properties of the 
time-dependent drift and diffusion processes with an absorbing barrier is pre- 
sented. For constant drift and diffusion, the conditional probabilities are usu- 



ally obtained via the method of images due to Lord Kelvin (see Feller, 1971 
p. 340). The applicability of this method to the solution of time-dependent 
problems and its limitations are discussed and a necessary and sufficient cri- 



terion is formulated in Section 2.2| The analysis is then extended to different 



functional forms of the time-dependent drift and diffusion terms. Section [3TT 
shows the analytical results for the first passage time statistics for a power-law 
time dependent drift and diffusion, while time-periodic drivers are analyzed 
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in Section 3.2 (see Jung 1993 Kim et al. 2010 Talkner et al. , 2005 and 



75 references therein, for a more comprehensive review of periodically-driven 

76 stochastic processes). Finally, in Section |4j we present a stochastic model of 

77 the total mountain water equivalent during the apex phase of the melting 

78 season, incorporating both temperature effects and snow-precipitation input 

79 in the form of a power-law time-dependent Bm with an absorbing boundary. 



so 2. Modeling Framework 

si When a time-dependent random forcing is the dominant driver of the dy- 

82 namics, a general representation for the state variable x(t) can be formulated 

83 in the form of a stochastic differential equation given by 

dx(t) = n(t) dt + a(t) dW{t) (1) 

84 where and a(t) are purely time-dependent drift and diffusion terms, 
as and W(t) is a Wiener process with independent and identically Gaussian 
se distributed (iid) increments W(t) - W(s) ~ 7V(0, t — s) for all t ^ s ^ 0. By 

87 assuming t = 0, the solution of Q takes the form 

x(t) = x + [ ^(s)ds+ [ [a(s)dW(s)\ (2) 
Jo Jo 

88 where t is time and Xq = x(0) can be either a random or a non-random 

89 initial condition independent of W{t) — 11/(0). The associated FP equation 

90 describing the evolution of the probability density function (pdf) of x{t) can 

91 be expressed as 

dp(x,t\x ) dp(x,t\x ) 1 2 d 2 p(x,t\x ) 
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92 where p(x,t\xo) is the transition pdf with initial condition 5(x — Xq) at to- 

93 Eq. ^ can also be expressed as a continuity equation for probability 

d . , . d ., . . fl . 

—p{x,t\x ) = -—j{x,t\x ), (4) 

94 where 

j(x,t|x ) = ii{t)p(x,t\x ) - -a (t) — , (5) 

95 is the probability current (or flux) and p(x, t\xo) is the conditional probability 

96 The solution of the FP equation in pi), is usually approached numerically 



97 (see for e.g., Schindler et al. (2005)). Whether this equation is analytically 

98 solvable for different functional forms of /i(t) and a(t) with an absorbing 

99 boundary, and whether these solutions can be applied in the study of the 

100 first passage statistics at such boundary is the main focus of this work. Case 

101 studies that employ these solutions are also presented. 

102 2.1. Solution with Natural Boundaries 

103 Consider first the solution of the FP equation ^ in the unbounded 

104 case. Given that the drift and diffusion coefficients depend only on time, 

105 the parabolic equation ^ can still be reduced to a constant-coefficient equa- 
loe tion of the form 

dp(z,r) d 2 p(z,r) 



dr dz 2 
107 by transforming the original variables x and t into 



(6) 



2 

and 



T = \ I a\t)dt + A (7) 



z = x- J fi(t)dt + B (8) 



109 where A and B are generic constants. The solution with natural boundaries 



is then (Polyanin 
condition 



2002) p(z,r) = exp ^— J^J. Hence, given the initial 



p(x,0\x ) = 6(x - x ), 



(9) 



the following normalized solution for an unrestricted process, starting from 
xq, can be obtained as 



p u (x,t\x ) 



1 



cxp 



(x-x - M(t)) 



21 



AS(t) 



(10) 



where, assuming the integrability of /i(t) and cr(t), 



and 



M(t) 



S(t) 



1 



fi(s)ds 



a 2 (s)ds. 



(11) 



(12) 



lie It should be noted that the transformation in equations and (m) also 

in applies to any boundary condition imposed at a finite position. Therefore, 

us as will be seen, it is not directly helpful in solving first passage time problems, 

n9 as in that case it would lead to a problem with moving absorbing boundary 

120 conditions. 

121 2.2. First Passage Time Distributions 

122 For a Bm process commencing at a generic position Xq at t — 0, the time 

123 at which this process reaches an arbitrary threshold a for the first time (first 

124 passage time) is itself a random variable whose statistics are fundamental in 

125 many branches of science such as chemistry, neural-sciences and economet- 

126 rics. In the following, it is assumed that the process is starting at a certain 
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state Xq > and that it is bounded to the positive semi-axis via an absorbing 
barrier x = 0. This hypothesis does not imply any loss of generality, consid- 
ering that the solution of Eq. ^ with an absorbing boundary condition only 
depends on the distance of the initial point Xq from the threshold, but not 
separately on x$ and the threshold position. Eq. (|3| is then solved with the 
boundary condition 

p(0,*) = 0, (13) 

and the additional condition of x = +oo being a natural boundary to ensure 
that j(+oo,t |xo) = 0. For such a system, the survival probability F(t \xq) 
is defined as the probability of the process trajectories not absorbed before 
time t, i.e. 



F(t \x ) 



p(x, t \xq) dx 



(14) 



and the first passage probability density g(t \x ) is either the "rate of de- 
crease" in time of F 



g(t\x ) 



d_ 
di 



F(t \x ) 



or, alternatively, the negative probability current at the boundary 

a 2 (t) 8 



g(t \x ] 



2 dx 



p(x,t\x ] 



\x=0 



(15) 



(16) 



wo since p(0 } t\x ) = from (13). 



mi 2. 3. Method of Images in Time-Dependent Systems 

142 When the drift and diffusion terms are independent of t and x, Eq. ^ 

143 with absorbing boundaries can be readily solved by the method of images, 



144 often adopted in problems of heat conduction and diffusion ( Cox & Miller 



us 1965 Daniels 1982; Lo et al. 2002; Redner 2001). This method can also 



146 be used for solving boundary- value problems for a Bm with particular forms 

147 of time- dependent drift and diffusion. The basic premise of this method is 

148 that given a linear PDE with a point source (or sink) subject to homogeneous 

149 boundary conditions in a finite domain, its general solution can be obtained as 

150 a superposition of many 'free space' solutions (i.e. disregarding the boundary 

151 conditions) for a number of virtual sources (i.e. outside the domain) selected 

152 so as to obtain the correct boundary condition. The image source (or sink) 

153 is placed as mirror image of the original source (or sink) from the boundary 

154 with a strength or intensity selected to match the boundary condition. 

155 Consider equation ^ with the conditions ^ and (13). To solve this 



156 problem with the method of images, the barrier at is replaced by a mirror 

157 source located at a generic point x = y, with y < such that the solutions of 

158 the Fokker-Planck equation emanating from the original and mirror sources 

159 exactly compensate each other at the position of the barrier at each instant 



wo of time (Redner, 2001). This implies the initial conditions in pj) must now 



161 be modified to 

p(x, 0) = S(x — xq) — exp (—77) S(x — y), (17) 

162 where 77 determines the strength of the mirror image source. Due to the lin- 

163 earity of the FP equation, the solution in the presence of the initial condition 



164 (17) is the superposition of elementary solutions 

p(x, t\x ) = p u {x, t\x ) - exp (-77) p u (x, t\y). (18) 



165 Since the condition (13) requires that p(0, t\xo) = 0, one obtains that 



(M(t)+x ) 2 _ , (M(t)+yf 

AS(t) ~ V+ AS(t) [ } 
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166 for all t > 0. By assuming t = 0, we have Xq = y 2 and recalling that y < 0, 



the resulting image position is — xq. This, inserted again in Eq. (19), yields 

M(t) 



S(t) 



rj_ 

X 



(20) 



168 where the constant q is analogous to the Peclet number of the process - i.e. 



169 the ratio between the advection and diffusion rates (Redner, 2001). 



170 After differentiating (20) with respect to t, it is seen that the method 



of images requires that the drift and the diffusion terms be proportional to 
each other. Namely, the intensity rj of the image source must be constant in 
time. In fact, only in this case it is still possible to transform the original 
time scale into a new one, for which the transformed process is governed 
by time-independent drift and diffusion terms. Hence, writing the drift and 
diffusion terms as 

'2 



jj,(t) = kh(t) and 



a- 



lh(t), 



177 the associated FP equation is 



dp 
di 



d d 2 
^ dx ^ ^ dx 2 



(21) 



(22) 



178 Transforming the original time t variable in 



T 



179 Equation (22) finally becomes 



dp 
df 



h(s)ds 



, 9 , d2 
-k— + l- 



(23) 



(24) 



dx dx 2 

This condition is valid for any time-dependent diffusion when the drift is 



i8i identically vanishing. Assuming the proportionality in (20) between fi(t) 



182 and cr(t), the general solution for ^ under conditions ([9J and (13) can be 

183 written as 



p{x,t\x ) 



cxp 



(x-XQ-M(t)) 2 



2^S(t) 

exp (— xoq) exp 



4S(i) 

(x+x -M(t)) 2 
4S(t) 



(25) 



184 provided M(t) = qS(t). Substituting for constant drift and diffusion in (25) 



185 one recovers the well-known solution for a biased Bm (|Cox & Miller| 1965) 

I J ™^ (xQ-X+flt) 



p(x,t\x ) 



2-K\/o 2 t 



exp 



2cr 2 t 



-exp (-^) exp 
we with survival function F(t|xo) given by 



(26) 



2x fi\ \ fit — x 



(27) 



o\ft j V cr 2 / [ 

187 where $ is the standard normal integral, and first passage time distribution 

(X + fitf 



g(t\xo) 



X 

aV2^tV 2 



exp 



2aH 



(28) 



las Equation (28) is the Wald (or inverse Gaussian) density function, that for a 

189 zero drift becomes of order t~ 3//2 as t — » +oo (the first passage time has no 

190 finite moments for pure diffusion). 

191 Similarly, the solution to the FP in equation (3) with a reflecting bound- 

192 ary at x = can be obtained by the method of images provided that drift 

193 and diffusion are proportional to each other. The solution then becomes 



P(x,t\x ) =^^{exp 



2^S(t) 

-!1f exp 



{x-x -M{t)) 2 
4S(t) 

1 — erf 



■ exp {—Xoq) exp 

\-x +M(t) \ 



(x+x -M(t)) 2 
4S(t) 



(29) 
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194 with | 



1 — erf 



x+x 



+M(t) \ 



2yJS(t) J 



being the Q-function representing the tail prob- 



195 ability of a Gaussian distribution. Equation (29) generalizes the solution in 



Cox & Miller (1965) for a Bm with constant drift and diffusion and a reflect- 



197 ing boundary at 0. 



198 3. Time Dependent Drift and Diffusion 

199 3.1. Power-Law Time Dependence 

200 As a first example of Bm with purely time-dependent drivers, the case of 

201 an unbiased diffusion (q = 0) and power-law time dependent diffusion term 

202 <J 2 (t) = 2At a and a > — 1 are considered. For this process, the conditional 

203 probability p(x,t\xo) with absorbing barriers at 0, takes on the form 

p(x,t\x ) = 

—exp 

204 while the survival function becomes 



VT±a f -(a+l) , 



t-(«+ 1 )(x-a:o) 2 (l+a) 
4A 



t -(q+i) ( x+Xo f : (l+a) 
4A 



(30) 



F{t\x c 



erf 



x (l + a) 2 
7= — 1 2 



(31) 



Figure [jja) shows the conditional probability ( p0| ) at a fixed time instance t = 
15 time steps for A = 15, xq = 50, and a = —0.1 (bold line) 0.5 (thin line), 
and 1 (dotted line). Given the asymptotic properties of the error function 



208 (Abramowitz & Stegun, 1964), the long-time behavior of F(t\x ) is then 



.to 



(1+q) 



2VA 



t 2 5 recovering for a = the —1/2 tail decay of an unbiased 



constant diffusion (see Figure [lib)). Also, by differentiating Eq. (31), one 
obtains 

In T Xn(a + l)r (a+1) " 



212 whose tail behaves as ~ t 



. Hence, Eq. (32) is an inverse Gaussian 



213 distribution - that for a = becomes an inverse Gamma distribution with 



214 shape parameter 1/2 (Johnson et al. 1994, pp. 284-285). These solutions 

215 characterize inter-arrival times between intermittent events when a system 



216 displays sporadic randomness (Gaspard & Wang, 1988; Molini et al. 2009 



Rigby & Porporato 2010). 



218 The solutions in the case of proportional power-law diffusion and drift can 

219 be derived in an analogous manner. For fj,(t) = qAt a and a(t) = s/2A l / 2 t a l 2 , 

220 the conditional probability p(x,t\xo) takes the form 



p(x,t\x ) = 
V5+Tt-("+ 1 )/ 2 



cxp 



4A 



— exp 



-qx - 



( 1+ a) t -C«+l)(x-^^+xo) : 



(33) 



-L4 



and the survival function, now incorporating the drift contribution, can be 
written as 



F(t\x ) 



(c+1) (Aqt a+1 +x +x a) 
l \ , T , f ■ (Aqt°'+ 1 ~x -xoa) 

—exp (— qx ) <P <t 2 ^ 



(34) 



2 X /A(a+1) 

For positive q's, F(t\xo) tends in the long term to 1 — exp(— qxo), while for 
negative g's, F(x,t\x ) ~ 2 ^^ 1 t~ SL ^ i exp ^~ 9 ^^px ~^ ■ This fact implies 
that the probability for a trajectory to be eventually absorbed is 1 for the 
biased process directed towards the barrier, and exp(— qxo) when the bias is 
away from the barrier (infinite aging). When the state variable represents 
the availability of a resource in time, the sign of q determines if this resource 
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is subject to continuos accumulation (positive q), or it undergoes a total 
depletion (negative q) with probability 1. Such a result is analogous to the 
one of a simple biased Bm with constant drift and diffusion (Redner 2001[ ), 
with the difference that in this case, F(t\x ) decays to or 1 — exp(— qx ) 
with a rate that is governed by a. 

As an example, Figures[l](c) and (d) respectively show a negatively biased 
power-law time-dependent Bm and a positively biased one for the same set 
of parameters in (b) and q = —0.1 and q = 0.1, for A = 1, xq = 1 and a = 
(constant diffusion, bold line), —0.5 (thin dotted line), 0.5 (dashed line), and 
1 (thin line). As evident in panel (c), F(x,t\x ) presents a faster decay to 
zero with increasing a, while for the positively biased Bm in panel (d) the 
decay to the asymptotic value 1 — exp(— gx ) is slower with decreasing a. 



24i Finally, g(t \xq) can be obtained from (34) as 



x l + a 3 / 2 

g(t\xo) = ^— 3+ Q exp 

2ViTAt^ 



j-(«+i) (Aqt a+1 +x + ax ) 



(35) 



4,4(1 + a) 

where for a = the decay of g(t\xo) recovers the constant diffusion t~ 3//2 -law 
for t — > oo and q = 0. 



244 3.2. Periodic Drift and Diffusion 

245 In this section, the case of a periodic diffusion in the form o~ 2 (t) = 

246 [2 A cos(cut)] 2 and q = is considered. For periodically driven diffusion, the 

247 conditional probability can be derived in the form 



p(x,t\x ] 



7T#(t) 

• exp 



exp 

ui(x-xo) 2 

m 



2w[x +xq) 



exp 



uj(x+xq) 2 

m 



(36) 
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where $(£) = A 2 [2ut + sm(2ut)] > 0. Thus, the solution becomes modulated 
in time with frequency u. The survival probability is in turn 



F(t\x ) = erf ( x 



to 



(37) 



that is represented in Figure [2]for different values of the frequency u. Finally, 
the first passage time density is an cu-modulated inverse Gaussian distribution 



. . . AxqA 2 u 3/2 cos(cut) 2 
g{t\x ) = — f= -Q/-tV , /9 — exp 



7T 



UXn 



In the case q ^ 0, the conditional probability p(x,t |x ) becomes 



(38) 



p(x,t\x ] 



— exp 



cxp 



-qx - 



u(xo+x- 



4oj , 



(39) 



where, again, the absorption at the barrier represents a recurrent (q < 0) or a 
transient (q > 0) state, as was observed for the power-law drift and diffusion 



255 process in Section 3.1 The recurrent case is illustrated in Figure (b) 



(d), where we report the time-position evolution of p(x, t\ function of 



increasing u. From (39), given j > 0, the expression for the survival function 



can be derived and takes the form 



F(t\x ) 



1 + erf 



q-d(t)+&x u 



+ exp (— qx ) erfc 



2 exp (-qx ) 



(40) 



which, given the equality erfc(— x) = 2 — erfc(x), can be alternatively ex- 
pressed as 



(41) 
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The first passage time density g{t\x ) is given by 



. . . 4A 2 x u) 3 2 cos(ut)'' 
g{t \x ) = — j= exp 



7T 



(q&(t) + Axpuf 
16w&(t) 



(42) 



The method of images can also be applied to the solution of different 
forms of periodic drivers, such as the case fj,(t) = q(B + Acos(cot)) and 
a(t) = \/2(B + Acos(ut)), with (B + Acos(ut)) > 0. In this last case, the 
drift term is the same as the one usually investigated in neuron dynamics 
by simple integrate-and-fire models displaying stochastic resonance (see for 



267 example the neuron dynamics case in Bulsara et al. 1994, 1995). In those 



models, the diffusion is usually constant so that the condition in equation 
(20) is not satisfied. Thus, it is often implied that fi(t) « a 2 / 2 to approxi- 
mately resemble a time dependent diffusion with drift identically vanishing or 
that B » A (approximating the simpler constant drift and diffusion case). 



Jn these cases, the method of images only offers approximated solutions (Bul- 



sara et al. , 11994, 1995)). Specifically, for a time dependent (and periodic) 



drift fi(t) = B + Acos(ut) and constant diffusion |cr , an approximation for 
p(x, t \xo ) in the presence of an absorbing barrier at can still be obtained by 
using the method of images conditional to the fact that fj,(t) « a 2 / 2. Only 
by adopting this assumption in fact, we can obtain an (approximated) solu- 



tion for the survival function by means of Eq. 25 although drift and diffusion 
are not strictly proportional to each other. In this way we find 

(43) 



F(t\x ) 



\ { erfc 



BU 



exp 



y/2cr<Jt 
2x (Buit+Asm(u)t)) 



a 2 ujt 



erfc 



Bt- 



A sin(ujt) 



-.X() 



and, analogous to Bulsara et al. (1994), from equation (15) the first passage 
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density can be expressed as 



xq exp 



[„. . Asin(ut) I 2 
\Bt+ ^--xq] 



A exp 



(44) 



282 The approximated nature of the solution is evidenced by the fact that, the 

283 image source intensity is no longer constant in time, so that by evaluating 

284 the probability current in we obtain 

[u(Bt - x ) + Asm(ut)} 2 



9(t\x Q ) 



Xq 
2W 3 / 2 



exp 



2oo 2 aH 



(45) 



which is different from (44). In any case, the first passage time pdf in equa- 



286 |f ion (44) is in good agreement with the numerical simulations in Bulsara 



et al. (1994 1995). Also, when A -> both the (44) and the (45) tend to 



288 the first passage time pdf for a simple biased Bm. 

289 As highlighted in Figure Q, when the magnitude of //(t) becomes sig- 

290 nificant, the two pdfs diverge due to the losses of probability density at the 



291 barrier (Eq. (45)). For this reason, the method of images cannot be con- 

292 sidered a general approach to solving problems described by Eq. ^ with a 

293 time-dependent Peclet number. 



294 4. A Case Study: Snowmelt Dynamics 

295 Snowmelt represents one of the paramount sources of freshwater for many 

296 regions of the world, and is sensitive to both temperature and precipitation 



297 ,fluctuations flBarnett et aLj |2004[ |2005[ |Barnett & Pierce[ |2009[ |Pepin k 



298 Lundquist, 2008; Perona et al. , 2001 You et al. 2010). Snow dynamics is 
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299 characterized by an accumulation phase during which snow water equivalent 

300 (i.e. the amount of liquid water potentially available by totally and instanta- 

301 neously melting the entire snowpack) increases until a seasonal maximum ho 

302 is reached, followed by a depletion phase in which the snow mantel gradually 

303 decays (and releases the stored water content) due to the increasing air tem- 

304 perature. Such a dynamics is complex and its general description requires 

305 numerous physical parameters that are rarely measured or available. In this 

306 section, we focus on a stochastic model describing the total water equivalent 

307 from both snow and rainfall during the melting season, as forced/fed by both 

308 precipitation (moving from the solid to the liquid precipitation regime) and 

309 increasing air temperature. 

310 Due to the simplified nature of our stochastic model, we will consider 
3n the total potential water availability (in terms of water equivalent) as the 
312 key variable, thus neglecting any further effects connected with snow per- 



313 eolation and metamorphism (De Walle & Rango, 2008). Snowfalls are here 

314 assumed to become more sporadic progressing into the warm season and the 

315 predominant controls over fresh water availability during the melting period 

316 are increasing air temperature and liquid precipitation. Accordingly, the 

317 melting phase is described by a power-law time dependent drift directed to- 

318 wards the total depletion of the snow mantle and by a power-law diffusion 

319 whose positive and negative excursions represent respectively precipitation 

320 events and pure melting periods. The melting process is often described by 

321 a linear function of time by using the so called "degree-day" approach with 



322 time- varying melting-rate coefficients (De Walle & Rango, 2008). Consid- 

323 ering that temperature varies seasonally and increases during the melting 
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324 season, a power-law form for drift and diffusion during the spring season, 

325 still represents a parametrically-parsimonious and effective approximation of 

326 the basic driver of the process. 

327 Under these assumptions, the dynamics of the total water equivalent 

328 depth for unit of area h - i.e. the amount of fresh water potentially available 



329 from both snow accumulation and rainfall (Bras, 1990) - at a given point in 

330 space can be can be reasonably described by the Langevin equation 

dh = -qkt a dt + V2ki^dW(t) (46) 

331 where k (with dimension L 2 /T a+1 ) represents the accumulation/ablation 

332 rate. Note that here h includes both the rainfall and snowmelt contribu- 

333 tions. Also, we hypothesize that both the drift and the diffusion scale with 

334 the same exponent a. This is a reasonable assumption given that variability 

335 of the process is expected to increase proceeding into the warm season. The 

336 initial condition is given by the snow water equivalent (SWE) ho, accumu- 

337 lated during the cold season. The survival probability F(t\ho) for a given 

338 initial SWE and the first passage time density g(t\h ) can be respectively 



339 calculated from rt34| and (35). Figure 5 shows few sample trajectories of 



340 the process (panel (a)) obtained by the numerical simulation of Eq. (46) 

341 by means of a forward Euler algorithm with a time step of 10 -2 days. The 

342 conditional probability p(h,t\h ) at different instants, the first passage time 

343 density g(t\h ), and the survival function F(t\h ), for the case a = 0.25 and 

344 k = 0.24 mm 2 /days a are also shown in panels (b) to (d). Here, we calibrated 

345 the parameters to obtain the mode of the first passage time at about 40 days 

346 after reaching the maximum SWE of the season h . The first passage time 

347 statistics presented offer important clues about the timing between melting 
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348 and summer fresh-water availability under different climatic scenarios (con- 

349 sider for example the FPT pdf in Figure (5^c)). 

350 5. Conclusions 

351 The first passage time properties of Brownian motion with purely time 

352 dependent drift and diffusion coefficients subjected to an absorbing barrier 

353 were investigated. These processes can be used to mimic a variety of en- 

354 vironmental and geophysical phenomena, representing the availability of a 

355 resource and its dynamics in time (e.g. the ablation phase of a snow mass 

356 accumulated during the winter period and forced by temperature and precip- 

357 itation). Survival functions and pdfs for the first passage times at the barrier 

358 were derived for power-law and periodic forcing time-dependent drift and 

359 diffusion terms for the associated Fokker Planck equation using the method 

360 of images. The general properties and limitations of this method were also 

361 reviewed, with reference to previous results obtained in the field of neural 

362 sciences and stochastic resonance. Particularly, we discussed how the ap- 

363 plicability of the method of images to a Bm with time-dependent drift and 

364 diffusion is limited to the case of a process with constant Peclet number, i.e. 

365 with a time-independent ratio of drift and diffusion. 

366 Where the time dependence is of the power-law type, the derived first 

367 passage time density and survival functions share many analogies with the 

368 statistics of inter arrival times between intermittent events when the con- 

369 sidered system displays sporadic randomness. In the case of a periodic 

370 time-dependence, first passage time statistics appear to be modulated by 

371 the frequency of the forcing. The periodic forcing case has been also used to 
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372 show the approximate nature of solutions obtained by the method of images, 

373 when time-dependent drift and diffusion terms are not linearly related. We 

374 finally show how a Bm with power-law decaying drift and diffusion can be 

375 used to describe the warm season dynamics of the total water equivalent in 

376 mountainous regions. 
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Figure 1: Conditional probability p(x,t\xo) at different fixed times t (a) and survival 
function F(t\xo) (b) for the pure power-law time dependent process described in Section 



3.1 together with F(t\xo) for the negatively biased power-law process (c) and for the 



positively biased one (d). Panel (a) represents p(x,t\xo) at a fixed time t — 15 steps for 
A = 15, .To = 50, and a = —0.1 (bold line), = 0.5 (thin line), and = 1 (dotted line). In (b) 
F(t\xo) is displayed as a function of t for A = 1, xo = 1 and a — (constant diffusion, bold 
line), a = —0.5 (thin dotted line), a — 0.5 (dashed line), and a = 1 (thin line). Panels 
(c) and (d) display respectively a negatively biased power-law time dependent Bm and a 
positively biased one for the same set of parameters in (b) and q = —0.1 and q = 0.1. 
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Figure 3: Conditional probability p(x, t\xo) for the periodic negatively biased Bm described 



in Section 3.2 Panel (a) represents p(x,t\xo) at a fixed time t = 3 steps for A = 15, 
x = 50, q = -0.05, and ui = 0.0001 (bold line), = 0.5 (thin line), and = 0.9 (dotted 
line). Also, contour plots (b) to (d) show p(x,t\xo) for A = 15, x^ — 450 and q = —0.01 
as a function of x and t, for lu = 0.0001 (panel (b)), u = 0.015 (panel (c)), and u> = 0.045 
(panel (d)) respectively. Note how the negative drift forces the probability mass toward 
the barrier. 
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Figure 4: First passage densities g(t\xo) (bold black line, Eq. 44 1, and g(t\xo) (red dotted 



line, Eq. 45), respectively for (a) (jt = 0.065, a = 0.5, x = 25, A = 0.032 and uj = 0.016; 



(b) (i = 0.065, a = 0.35, x Q = 15.5, A = 0.025 and w = 0.04; (c) /i = 0.065, cr = 0.2, 
.t = 25, A = 0.03 and w = 0.07, and (d) n = 0.065, cr = 0.2, x n = 25, A = 0.03 and 
w = 0.15. The discrepancy between g(t\xo) and g(t\%o) clearly signifies the failure of the 
method of images for problems with time-dependent Peclet numbers. 
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Figure 5: Sample trajectories of specific water equivalent from elevated regions during the 
melting season (a), and analytical results for the coupled stochastic melting-precipitation 



process of equation (46) (panels (b) to (d)). Numerical results were obtained by simulating 
Eq. (46) by means of an Euler algorithm with step 10 -2 days. Panel (a) shows few 
sample trajectories of the process together with the curve of maximum values (upper curve) 
and minimum values (lower curve) over an ensemble of 10000 simulations, for a = 0.25 
and k = 0.24 mm 2 /days Q . Analytical results for the conditional probability p(h,t\ho) at 
different instants, the first passage time density g(t\ho), and the survival function F(t\ho), 
are also shown in panels (b) to (d). 
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